Eigenfunction Expansion Based Galerkin Approaches for Stochastic Finite Element Analysis

A novel approach is suggested to compute the response of discretized stochastic elliptic partial differential equations by utilising the stochastic finite element analysis method. The mathematical form of the approach is established by projecting random scalars onto a random basis. By implementing an eigendecomposition of a system’s stiffness matrix, the random scalars and random basis are computed. Computational reduction is achieved by approximating the random eigensolutions and by only including dominant terms. Two novel error minimisation techniques have been proposed in order to lower the error introduced by the approximations and the truncation: i) A weak Galerkin approach, ii) A strong Galerkin approach. The proposed methods are applied to analyse the bending of a stochastic cantilever beam. The results obtained through the proposed approaches are compared with those obtained by using direct Monte Carlo Simulations and by using polynomial chaos. 1 INTRODUCTION Uncertainties can substantially affect the analysis of physical structures. These uncertainties can occur in the properties of the material, in the geometry or boundary conditions of the structure or in the applied loads [1]. In order to represent the uncertainties that occur in physical systems, a stochastic finite element method [SFEM] can be applied. In this work, a stochastic elliptic partial differential equation is considered − ∇ , ∇ , = (1) with the associated Dirichlet condition , = 0 (2) The spatial dimension under consideration is a bounded domain ∈ R with piecewise Lipschitz boundary , where is less than four. Ω, F, ) is a probability space where ∈ Ω is a sample point from the sampling space Ω, F is the complete -algebra over the subsets of Ω and is the probability measure. In Equation (1) ∶ R × Ω → R is a random field [2], which can be viewed as a set of random variables indexed by ∈ R . We assume the random field , to be stationary and square integrable. , is also able to model different physical quantities. Through combining SFEM with a Karhunen-Loève expansion, the discretized equivalent of Equation (1) can be been obtained. The details of obtaining the S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 75 discretized equivalent of Equation (1) has been omitted, but can be located in numerous textbooks including [3]. For the case of a static cantilever beam experiencing a deterministic applied force, the discretized equivalent can take the following form !" + $%& !& ' &() * = ! = + (3) where !" ∈ R,×, represents a deterministic, positive definite, symmetric matrix. !∈ R,×, are random symmetric matrices for i = 1, 2, ... ,M, ∈ R, the response vector, ξ& a set of random variables and + ∈ R, the deterministic input force vector. For the remainder of this paper we will assume ξ& to be Gaussian random variables. We aim to produce a new solution approach to solve the discretised equivalent of Equation (1). 2 MOVIVATION The exact solution to the set of stochastic linear equations given above can be obtained through a direct Monte Carlo Simulations [DMCS] approach by solving the following expression for each sample. = ! 45 + (4) However, DMCS can be seen as a computationally expensive method, especially if there is a large number of stochastic linear equations to be solved. In order to avoid the use of DMCS, alternative methods have been explored. The response of Equation (3) can be represented through summing products of random scalars and deterministic vectors =$ 6 76 '8 6() (5) Equation (5) can be considered as the polynomial chaos method [PC] where 6 would correspond to the polynomial chaoses, and 76 would correspond to unknown deterministic vectors. The full details of this method can be found in [3]. Due to the vector in Equations (5) being deterministic, we aim to acquire a solution where both the scalars and vectors are random. =$ 6 76 '9 6() (6) We aim to see if Equation (6) can incorporate the full stochastic nature of Equation (3). Therefore the aim of this paper is to obtain an expression for the response of Equation (3) that is of the same form as Equation (6). 3 INTRODUCING THE METHOD In order to implement our aim, the random eigenvalue problem is initially considered ! :; = <; :; ; > = 1, 2, ... (7) For convenience, the matrices of the random eigenvalues and eigenvectors of ! are defined as follows S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 76 ? = diag <) , < ,... , <, ∈ R,×, (8) C = :) ,: ,... ,:, ∈ R,×, (9) The random eigenvalues are arranged in ascending order so <) < < < ⋯ < <, . The corresponding eigenvectors are arranged in the same order. Due to the orthogonality of C , we can deduce that C 45 = C F. Thus the following identities can be defined ( has been omitted for notational convenience) CG! C = ?; ! = C4F? C45 and !4) = C ?45 CG (10) Using these identities, the response of Equation (3) can be expressed as = C ?45 CG + = $:6G + <6 , 6() :6 (11) Equation (11) is of the same form as Equation (6) where :IJ K + LI K corresponds to the scalar term 6 , :6 corresponds to the vector term 76 and corresponds to M . In this particular method, :6 forms a complete orthogonal basis. Therefore, it can be concluded that the response of Equation (3) can be expressed in the same form as Equation (6) where we have random scalars projected onto a stochastic basis. For the remainder of this paper, this method has been labelled [SP]. 4 COMPUTATIONAL REDUCTION 4.1 Approximating the random eigenvalues and eigenvectors Approximating the random eigenvalues and eigenvectors may improve the calculation cost and there are numerous methods of doing so. Due to its low computational cost, a first order perturbation approach for obtaining the random eigenvalues and eigenvectors has been explored. Solutions of different perturbation methods are obtained by truncating the Taylor series expansion. Due to its efficiency and ease, the first order perturbation method has been used. The Nth random eigenvalue and its corresponding random eigenvector is given by <6 ≈ <6P + ∑ RSLI STUV %; ';() and :6 ≈ :6P + ∑ RS:I STUV %; ';() (12) where %; is a set of Gaussian random variables with mean zero and unit variance. By differentiating the eigenvalue equation with respect to %;, pre-multiplying with :6G and utilising that :6G :6 = 1, the derivatives of the eigenvalues and eigenvectors can be expressed as <6 %; = :6P G ! %; :6P and :6 %; = $ :&P; W ! %; :6P <6P − <&P; :; X &() Y6 (13) where S! STU = !;. The derivative of the eigenvectors requires all the deterministic eigenvalues and eigenvectors to be known. However care is needed with this approach as all the eigenvalues need to be unique and the coefficient of variation needs to be of a moderate S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 77 value. For the case of repeated eigenvalues, the general method of approximating the response of Equation (3) continues to be valid. However a different approach would be needed whilst approximating the eigensolutions to that given in this section. 4.2 Truncation The series given in Equation (11) can be truncated after a certain amount of terms. The high terms of the summation have a relatively low value due to the eigenvalues being ordered ascendingly; this allows the low valued terms to be discarded whilst retaining the dominant terms in the series. Consequently Equation (11) can be truncated as follows ≈$:6G + <6 Z 6() :6 (14) where <6 and :6 represent the random eigenvalues and the random eigenvectors and t corresponds to the number of terms retained in the summation. A low-cost MCS is performed in order to obtain the full response of . 5 ERROR MINIMISATION The approximations and truncation introduced in Section 4 introduces error into the calculation. This has motivated an error minimisation appraoch, and as a consequence, two Galerkin approaches are considered: (a) a weak Galerkin approach [SPWG] (b) a strong Galerkin approach [SPSG]. The SPWG approach is initially considered. 5.1 A weak Galerkin approach [SPWG] For this approach, the solution vector is modified to take the following form [ ≈$\6 ]:6G + <6 ^ Z 6() :6 (15) where <6 ∈ R and :6 ∈ RX×X represent the approximated random eigenvalues and eigenvectors, + ∈ RX×X the deterministic applied force and \6 ∈ RZ are deterministic constants which need to be determined. The residual vector for this the new approach is defined as _) = ! [ − + (16) By making the residual orthogonal to a basis function, the deterministic scalars \6 can be computed. As Equation (15) can be viewed as a projection onto a subset of the random eigenvectors, the residual can be made orthogonal to the same subset of random eigenvectors. _) ⊥ :; ∀ > = 1, 2, ... , b (17) Here 〈d ,e 〉 = g d e h K defines the inner product norm. By using this condition and the expression for the residual, one has i j:;W kl ω n$\6 ]:6G + <6 ^ Z 6() :6 o − +pq = 0 (18) S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 78 where ir∎t donates the expected value. For notational convenience, we can define u6 = v:IJ K + LI K w, thus it can be shown that Equation (18) can take the following form ix$:;W l :6 u6 \6 Z 6() y = ir:; t+ (19) By defining the vector z = \), \ , ... , \Z G, Equation (19) can be simplified to ir{ t z = ir| t N, > = 1, 2, ... , b (20) where {;6 = :;W l :6 u6 ; ∀ N, > = 1,2, ... , b, and | = :; +; ∀ > = 1,2, ... , b. Therefore by solving the set of linear equations given by Equation (20), an explicit closed form for the unknown coefficients can be obtained. The arising expected values can be computed by using low-order fast Monte Carlo Simulations. 5.2 A strong Galerkin approach [SPSG] In a similar manner to the previous approach, the solution vector can be modified to take the following form ́ ≈$ 6 ]:6G + <6 ^ Z 6() :6 (21) where <6 ∈ R and :6 ∈ RX×X represent the approximated random eigenvalues and eigenvectors, ~ ∈ RX×X the deterministic applied force. Contrary to the previous approach, 6 ∈ R are unknown constants that need to computed for each realisation. The residual vector for this approach is defined as _ = ! ́ − + (22) By making the residual orthogonal to a basis function, 6 can be computed. By using the same analogy as seen in the previous approach, the residual can be made orthogonal to the subset of random eigenvectors. This results in the following expression :;W kl ω n$ 6 ]:6G + <6 ^ Z 6() :6 o − +p = 0 (23) For notational convenience, we can define u6 = v:IJ K + LI K w. Therefore, the above expression can be manipulated to give $:;W l :6 u6 6 Z 6() € = :; + (24) where 6 would be computed for each realisation. By defining the vector  = ) , ,... , Z G, Equation (15) can be simplified to {  = | N, > = 1, 2, ... , b (25) S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 79 where {;6 = :;W l :6 u6 ; ∀ N, > = 1,2, ... , b, and | = :; +; ∀ > = 1,2, ... , b. The number of equations that need to be solved in order to calculate the unknown vector  corresponds to the value of b. Therefore, similarly to the weak Galerkin approach, the lower the dimension of the reduced system, the fewer the number of equations that need to be solved. 6 APPLICATION Thus far we have discussed five different methods for computing or approximating the response of Equation (3): i) DMCS, ii) SP, iii) SPWG, iv) SPSG and v) PC. These five methods are now applied to compute the bending of a static Euler-Bernoulli cantilever beam of length 1.00 m. The cross-section of the beam is a rectangle of length 0.04 m and height 0.005 m and a 1.00 N deterministic vertical point load is applied at the free end of the beam. The system has been discretized into a 100 elements by using SFEM. Consequently, the dimension of the corresponding determinant matrix is 200 × 200. For the deterministic case, the Young's modulus is E = 200 × 10 Nm thus corresponding to a steel beam. The deterministic second moment of area (moment of inertia) of the beam is 4.1667 × 10 m. The bending rigidity, ‚ƒ, can be assumed to be a stationary Gaussian random field ‚ƒ „, = ‚ƒ ……… †1 + „, ‡ (26) where ‚ƒ ……… = 83.33 Nm. The function „, represents a stationary Gaussian field with zero mean, with „ being the coordinate direction of the length of the beam. The standard deviation is given by 0.18 × ‚ƒ ………, and the covariance function by ‹Œ „), „ = Œ  |849| / ‘’ (27) where “Œ is the correlation length (0.50 m) and Œthe standard deviation; The KL expansion of the system has been truncated to include two terms. The five listed methods have been simulated 10,000 times. For the SP, SPWG and SPSG methods, 4 terms have been retained in their respective equations, and a 4 order polynomial chaos approach is used for the PC method. Figures 1a and 1b illustrate the mean and standard deviation of the normalised vertical displacement at all nodes of the beam. It is apparent that all the approximation methods captures the mean of the DMCS method pretty well. The SPSG method best captures the standard deviation of the DMCS method, whilst both the SP and SPWG methods considerably out performs the PC approach. The percentage error of the vertical displacement arising when using the SP, SPWG, SPSG and the PC methods in place of the DMCS method is illustrated in Figures 1c and 1d. The percentage error is given by ”% = 100 × |–M‹— − ‹M| –M‹— (28) where DMCS indicates the solution of the direct approach, and CM the solution of the comparable methods. Barring the initial 0.20m of the bar, the percentage error of the mean and standard deviation of the SP, SPWG and SPSG methods are considerably lower than that of the PC method. When comparing the SPWG and SPSG methods it is apparent that the percentage error of the mean is always slightly lower when the SPSG method is used. However when comparing the percentage error of the standard deviation, barring the initial 0.30m of the bar, the SPSG method considerably outperforms the SPWG method. S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 80 (a) Mean of the vertical displacement (b) Standard deviation of the vertical displacement (c) Percentage error of the mean (d) Percentage error of the standard deviation The ˜ relative error for the mean and standard deviation of the response of the cantilever beam is considered. The ˜ relative error of the mean is defined as ”™9 = ‖›œ'ž − ›' ‖ ‖›œ'ž‖ (29) where ›œ'ž denotes the mean of the response vector obtained by using the DMCS method and ›' the mean of the response vector of the comparable methods. The expression for the ˜ relative error of the standard deviation takes a similar form to that of Equation (29). Table 1 displays the ˜ relative error for the mean and the standard deviation for different truncation values of Equations (11), (15) and (21) and for the PC approach. It is apparent that both the SPWG and SPSG methods produces lower relative errors than the SP method for the mean. When considering the ˜ relative error for the standard deviation, the SPWG method does not significantly lower the error. However a substantial drop in the relative error is seen when the SPSG method is used. As the relative errors for both the mean and standard deviation decreases as the number of terms in the summations increase, it is apparent that the truncation value can be chosen in line with an acceptable level of error. The CPU times for the five methods is given in Table 2. Table 1: ˜ error arising in the mean and standard deviation of the response for each of the methods Relative error in the mean Relative error in the standard deviation Number of terms SP SPWG SPSG PC SP SPWG SPSG PC 2 0.0063 0.0058 0.0058 0.0092 0.0092 0.0029 4 0.0036 9.59 × 10 9.07 × 10 0.0101 0.0101 0.0010 5 0.0025 6.37 × 10 5.54 × 10 0.0102 0.0101 5.14 × 10-4 10 0.0025 3.28 × 10 9.54 × 10 0.0102 0.0102 9.88 × 10-5 0.0017 0.0635 S. E. Pryse and S. Adhikari Meccanica dei Materiali e delle Strutture | VI (2016), 1, PP. 74-81 81 Table 2: Computational time Method DMCS SP SPWG SPSG PC CPU (sec) 10.21 0.76 3.25 3.24 0.84 Regarding CPU, the SP method outperforms all the other methods. Although that using the SPWG and SPSG methods is faster than the DMCS method, the PC method is faster. However when the order of the PC method or the number of degrees of freedom is sufficiently increased, the CPU of the PC method will substantially increase and in turn, will be higher than the SPWG and SPSG methods. Generally it can be concluded that the SPSG method outperforms the SPWG method as it produces less error in a similar CPU time. However the SPWG method does not need as much storage capacity as the SPSG method. 7 CONCLUSION An approach has been suggested to calculate the response of discretized stochastic elliptic partial differential equations. Through utilising the stochastic finite element method and the random eigenvalue problem, it has been proven that the solution can be represented through projecting random scalar onto a random basis. Due to the high computational cost associated with calculating the exact solution, a reduced approach is proposed where random eigenvalues and eigenvectors are approximated and low valued terms are discarded. Two novel multiplicative Galerkin error minimisation approaches have been presented. i) a weak Galerkin approach, ii) a strong Galerkin approach. The proposed methods have been used to analyse the bending of a static cantilever beam. The solutions obtained through the proposed methods have been compared with those obtained through direct MCS and through using a polynomial chaos method. Although both error minimisation approaches lowers the error of the mean, it is only the strong Galerkin error minimisation approach that substantially lowers the error arising in the standard deviation. Further work would focus on efficient ways of computing random eigenfunctions and performing model-order reduction. 8 ACKNOWLEDGEMENTSThe authors acknowledge the financial support received from Engineering Research NetworkWales (one of three Sêr Cymru National Research Networks). BIBLIOGRAPHY[1] G. Maymon, Structural Dynamics and Probabilistic Analysis for Engineers,Butterworth-Heinemann, (2008).[2] E. VanMarcke, Random Fields: Analysis and Synthesis, MIT Press, Cambridge MA,(1983).[3] R. G. Ghanem and P. D. Spanon, Stochastic Finite Elements: A Spectral Approach(revised edition), Dover Publications Inc., (2012).